####Implicit Attitudes Towards an Authoritarian Regime - JOP - R Analysis File####
###Rory Truex, Assistant Professor, Princeton University Department of Politics and Woodrow Wilson School of Public and International Affairs, rtruex@princeton.edu###
###Daniel Tavana, Ph.D Candidate, Princeton University Department of Politics

setwd('/Users/rtruex/Google Drive/Implicit Attitudes Towards Authoritarian Regime/')
rm(list=ls(all=TRUE))

##LOAD DATA AND PACKAGES##

install.packages("arm")
install.packages("foreign")
install.packages("Zelig")
install.packages("ggplot2")
install.packages("gplots")
install.packages("RColorBrewer")
install.packages("Hmisc")
install.packages("list")
install.packages("plyr")
install.packages("ebal")
install.packages("doBy")

library("arm")
library("foreign")
library("Zelig")
library("ggplot2")
library("gplots")
library("RColorBrewer", lib.loc="/Library/Frameworks/R.framework/Versions/3.1/Resources/library")
library("Hmisc")
library("list")
library("plyr")
library("ebal") 
library("doBy")

require(grid)
vp.layout <- function(x, y) viewport(layout.pos.row=x, layout.pos.col=y)
arrange_ggplot2 <- function(..., nrow=NULL, ncol=NULL, as.table=FALSE) {
  dots <- list(...)
  n <- length(dots)
  if(is.null(nrow) & is.null(ncol)) { nrow = floor(n/2) ; ncol = ceiling(n/nrow)}
  if(is.null(nrow)) { nrow = ceiling(n/ncol)}
  if(is.null(ncol)) { ncol = ceiling(n/nrow)}
  
  grid.newpage()
  pushViewport(viewport(layout=grid.layout(nrow,ncol) ) )
  ii.p <- 1
  for(ii.row in seq(1, nrow)){
    ii.table.row <- ii.row  
    if(as.table) {ii.table.row <- nrow - ii.table.row + 1}
    for(ii.col in seq(1, ncol)){
      ii.table <- ii.p
      if(ii.p > n) break
      print(dots[[ii.table]], vp=vp.layout(ii.table.row, ii.col))
      ii.p <- ii.p + 1
    }
  }
}

data<-read.csv("data.final.csv") 
iat<-read.csv("iat.csv")
list2a<-read.csv("list2a.csv") 

##DATA CLEANING## 

data<-merge(data,list2a,by="session_id",all=TRUE)
data[data=="98"]<-NA
data<-subset(data, data$filter.complete==0 & data$filter.error==0 & data$filter.calculated==0)

data$female<-data$gender-1

data$birthyear[data$birthyear=="99"]<-NA
data$age <- data$birthyear+15

for(level in unique(data$province)){
  data[paste("province.ind", level, sep = "")] <- ifelse(data$province == level, 1, 0)
}

data$location <- data$province
data$location[data$location=="1"]<-"alexandria"
data$location[data$location=="2"]<-"aswan"
data$location[data$location=="3"]<-"asyut"
data$location[data$location=="4"]<-"beheira"
data$location[data$location=="5"]<-"benisuef"
data$location[data$location=="6"]<-"cairo"
data$location[data$location=="7"]<-"dakahlia"
data$location[data$location=="8"]<-"damietta"
data$location[data$location=="9"]<-"fayoum"
data$location[data$location=="10"]<-"gharbia"
data$location[data$location=="11"]<-"giza"
data$location[data$location=="12"]<-"ismailia"
data$location[data$location=="13"]<-"kafrelsheikh"
data$location[data$location=="14"]<-"luxor"
data$location[data$location=="15"]<-"matruh"
data$location[data$location=="16"]<-"minya"
data$location[data$location=="17"]<-"monufia"
data$location[data$location=="18"]<-"newvalley"
data$location[data$location=="19"]<-"northsinai"
data$location[data$location=="20"]<-"portsaid"
data$location[data$location=="21"]<-"qalyubia"
data$location[data$location=="22"]<-"qena"
data$location[data$location=="23"]<-"redsea"
data$location[data$location=="24"]<-"sharqia"
data$location[data$location=="25"]<-"sohag"
data$location[data$location=="26"]<-"southsinai"
data$location[data$location=="27"]<-"suez"
data$location[data$location=="98"]<-NA

data$region <- data$location
data$region[data$region=="alexandria"]<-"alex"
data$region[data$region=="aswan"]<-"upper"
data$region[data$region=="asyut"]<-"upper"
data$region[data$region=="beheira"]<-"delta"
data$region[data$region=="benisuef"]<-"upper"
data$region[data$region=="cairo"]<-"cairo"
data$region[data$region=="dakahlia"]<-"delta"
data$region[data$region=="damietta"]<-"delta"
data$region[data$region=="fayoum"]<-"upper"
data$region[data$region=="gharbia"]<-"delta"
data$region[data$region=="giza"]<-"cairo"
data$region[data$region=="ismailia"]<-"delta"
data$region[data$region=="kafrelsheikh"]<-"delta"
data$region[data$region=="luxor"]<-"upper"
data$region[data$region=="matruh"]<-"desert"
data$region[data$region=="minya"]<-"upper"
data$region[data$region=="monufia"]<-"delta"
data$region[data$region=="newvalley"]<-"desert"
data$region[data$region=="northsinai"]<-"desert"
data$region[data$region=="portsaid"]<-"delta"
data$region[data$region=="qalyubia"]<-"delta"
data$region[data$region=="qena"]<-"upper"
data$region[data$region=="redsea"]<-"desert"
data$region[data$region=="sharqia"]<-"delta"
data$region[data$region=="sohag"]<-"upper"
data$region[data$region=="southsinai"]<-"desert"
data$region[data$region=="suez"]<-"delta"

for(level in unique(data$region)){
  data[paste("region.ind", level, sep = "")] <- ifelse(data$region == level, 1, 0)
}

data$lowed<-NA
data$lowed[data$edu>4]<-0
data$lowed[data$edu<5]<-1

data$charitable[data$charitable==2]<-0

data$prof_assoc[data$prof_assoc==2]<-0

data$work_employed<-NA
data$work_employed[data$work<3]<-1
data$work_employed[data$work>2]<-0

data$work_student<-0
data$work_student[data$work==5]<-1
data$work_student[data$work==NA]<-NA

data$work_retired<-0
data$work_retired[data$work==3]<-1
data$work_retired[data$work==NA]<-NA

data$work_govemp<-0
data$work_govemp[data$work_sector==1]<-1

data$rel_christian<-0
data$rel_christian[data$religion==2]<-1
data$rel_christian[data$religion==NA]<-NA

data$rel_muslim<-0
data$rel_muslim[data$religion==1]<-1
data$rel_muslim[data$religion==NA]<-NA

data$expenses<-data$spend_monthly

data$highinc<-0
data$highinc[data$expenses>6]<-1
data$highinc[is.na(data$expenses)==TRUE]<-NA

data$vote_parliament[data$vote_parliament==2]<-0

data$vote_president[data$vote_president==2]<-0

data$approve_abdel_rev<-5-data$approve_abdel

data$direct.support<-0
data$direct.support[data$approve_abdel_rev==4]<-1
data$direct.support[data$approve_abdel_rev==3]<-1
data$direct.support[is.na(data$approve_abdel_rev)==TRUE]<-NA

data$implicit.support<-0
data$implicit.support[data$STIAT.test>0]<-1
data$implicit.support[is.na(data$STIAT.test)==TRUE]<-NA
mean(data$implicit.support)

data$implicit.support.15<-0
data$implicit.support.15[data$STIAT.test>.15]<-1
data$implicit.support.15[is.na(data$STIAT.test)==TRUE]<-NA
mean(data$implicit.support.15)

data$implicit.against.15<-0
data$implicit.against.15[-0.15>data$STIAT.test]<-1
data$implicit.against.15[is.na(data$STIAT.test)==TRUE]<-NA
mean(data$implicit.against.15)

data$poldemoc <- rowMeans(subset(data, select = c(democracy1, democracy3)), na.rm = TRUE)
data$poldemoc <- 5-data$poldemoc

data$polislam <- rowMeans(subset(data, select = c(democracy4, democracy5)), na.rm = TRUE)
data$polislam <- 5-data$polislam

data$islamist<-0
data$islamist[data$polislam >2.5]<-1
data$islamist[data$religion==NA]<-NA
data$islamist[data$rel_muslim==0]<-0
data$islamist[is.na(data$polislam)==TRUE]<-NA
sum(data$islamist, na.rm=TRUE)

data$liberal<-0
data$liberal[data$poldemoc>2.5 & data$polislam<2]<-1
data$liberal[data$islamist==1]<-0
data$liberal[is.na(data$polislam)==TRUE]<-NA
data$liberal[is.na(data$poldemoc)==TRUE]<-NA
sum(data$liberal, na.rm=TRUE)

##ANALYSIS##

#Table A1: Sample Comparisons

#Implement Weighting Using Entropy Balancing

ab.data <- read.csv("abw3.csv", header = TRUE) 
ab.data <- ab.data[ab.data$country == "Egypt", ]
ab.data <- ab.data[ab.data$q1011a == "Yes", ]

ab.data$female<-0
ab.data$female[ab.data$sex=="Female"]<-1

ab.data$age<-ab.data$q1001

ab.data$rel_christian<-0
ab.data$rel_christian[ab.data$q1012=="Christian"]<-1
ab.data$rel_christian[ab.data$q1012=="Missing"]<-NA

ab.data$lowed<-1
ab.data$lowed[ab.data$q1003=="ba"]<-0
ab.data$lowed[ab.data$q1003=="MA and above "]<-0

ab.data$work_employed<-0
ab.data$work_employed[ab.data$q1004=="Yes"]<-1

ab.data$work_student<-0
ab.data$work_student[ab.data$q1005=="A student"]<-1

ab.data$work_retired<-0
ab.data$work_retired[ab.data$q1005=="Retired"]<-1

ab.data$work_govemp<-0
ab.data$work_govemp[ab.data$q1007=="A governmental employee"]<-1
ab.data$work_govemp[ab.data$q1007=="Director of an institution or a high ranking governmental employee"]<-1

ab.data$charitable<-0
ab.data$charitable[ab.data$q5012=="Yes"]<-1

ab.data$prof_assoc<-0
ab.data$prof_assoc[ab.data$q5013=="Yes"]<-1

ab.data$highinc<-0
ab.data$highinc[ab.data$q1015>3000]<-1
ab.data$highinc[ab.data$q1015>50000]<-NA

provinces<-unique(ab.data$q1)

ab.data$region<-NA
ab.data$region[ab.data$q1=="The West Bank"]<-"cairo"
ab.data$region[ab.data$q1=="Cairo"]<-"cairo"
ab.data$region[ab.data$q1=="Suez"]<-"delta"
ab.data$region[ab.data$q1=="Damietta"]<-"delta"
ab.data$region[ab.data$q1=="Dakahlia"]<-"delta"
ab.data$region[ab.data$q1=="El Sharqiya"]<-"delta"
ab.data$region[ab.data$q1=="Qalyubia"]<-"delta"
ab.data$region[ab.data$q1=="Monufia"]<-"delta"
ab.data$region[ab.data$q1=="Giza"]<-"cairo"
ab.data$region[ab.data$q1=="Faiyum"]<-"upper"
ab.data$region[ab.data$q1=="Qena"]<-"upper"
ab.data$region[ab.data$q1=="Alexander"]<-"alex"
ab.data$region[ab.data$q1=="Minya"]<-"upper"
ab.data$region[ab.data$q1=="The Red Sea"]<-"desert"
ab.data$region[ab.data$q1=="Kafir el-Sheikh"]<-"delta"
ab.data$region[ab.data$q1=="Sohag"]<-"upper"
ab.data$region[ab.data$q1=="Beheira"]<-"delta"
ab.data$region[ab.data$q1=="Ismailia"]<-"delta"
ab.data$region[ab.data$q1=="Asyut"]<-"upper"
ab.data$region[ab.data$q1=="Luxor"]<-"upper"
ab.data$region[ab.data$q1=="Beni Suef"]<-"upper"
ab.data$region[ab.data$q1=="Aswan"]<-"upper"

for(level in unique(ab.data$region)){
  ab.data[paste("region.ind", level, sep = "")] <- ifelse(ab.data$region == level, 1, 0)
}

#Implement Entropy Balancing
ab.data.bal<-data.frame(cbind(ab.data$bid, ab.data$female, ab.data$age, ab.data$rel_christian, ab.data$lowed, ab.data$work_employed, ab.data$work_student, ab.data$work_retired, ab.data$work_govemp, ab.data$charitable, ab.data$prof_assoc, ab.data$region.inddelta, ab.data$region.indupper, ab.data$region.indalex, ab.data$region.inddesert))
colnames(ab.data.bal)<-c("session_id","female","age","rel_christian","lowed","work_employed","work_student","work_retired","work_govemp","charitable","prof_assoc","region.inddelta", "region.indupper", "region.indalex", "region.inddesert")
ab.data.bal$treat<-1
sciat.data.bal<-data.frame(cbind(data$session_id, data$female, data$age, data$rel_christian, data$lowed, data$work_employed, data$work_student, data$work_retired, data$work_govemp, data$charitable, data$prof_assoc, data$region.inddelta, data$region.indupper, data$region.indalex, data$region.inddesert))
colnames(sciat.data.bal)<-c("session_id","female","age","rel_christian","lowed","work_employed","work_student","work_retired","work_govemp","charitable","prof_assoc","region.inddelta", "region.indupper", "region.indalex", "region.inddesert")
sciat.data.bal$treat<-0

data.bal<-rbind(sciat.data.bal, ab.data.bal)
data.bal<-na.omit(data.bal)
treat<-data.bal$treat

X<-cbind(data.bal$female,data.bal$age, data.bal$rel_christian, data.bal$lowed, data.bal$work_employed, data.bal$work_student, data.bal$work_retired, data.bal$work_govemp, data.bal$charitable, data.bal$prof_assoc, data.bal$region.inddelta, data.bal$region.indupper, data.bal$region.indalex, data.bal$region.inddesert)
eb.out<-ebalance(Treatment=treat, X=X, base.weight = NULL, norm.constant = NULL, coefs = NULL, max.iterations = 500, constraint.tolerance = 1, print.level = 0)

apply(X[treat==1,],2,mean)
apply(X[treat==0,],2,weighted.mean,w=eb.out$w)
apply(X[treat==0,],2,mean)

eb.out.tr <- ebalance.trim(eb.out, print.level=3, max.weight=20)
apply(X[treat==0,],2,weighted.mean,w=eb.out.tr$w)

weights<-subset(data.bal, data.bal$treat==0)
weights$w<-eb.out$w
weights$w.tr<-eb.out.tr$w
weights<-cbind(weights$session_id, weights$w,weights$w.tr)
colnames(weights)<-c("session_id","w","w.tr") 

data<-merge(data,weights,by="session_id",all=TRUE)
data$w[is.na(data$w)==TRUE]<-0

hist(data$w.tr)

#Figure 3: Distribution of D-score

pdf('fig1-dscore-hist.pdf', width=6, height=4)                                                                                                                                                                                                                                                                                                                                                                                                                                              
ggplot() + geom_rect(data=NULL,aes(xmin=0,xmax=.15,ymin=-Inf,ymax=Inf),fill="darkred",  alpha=0.1) + geom_rect(data=NULL,aes(xmin=0.65,xmax=Inf,ymin=-Inf,ymax=Inf),fill="darkred",  alpha=0.4) + geom_rect(data=NULL,aes(xmin=0.35,xmax=.65,ymin=-Inf,ymax=Inf),fill="darkred", alpha=0.3) + geom_rect(data=NULL,aes(xmin=0.15,xmax=.35,ymin=-Inf,ymax=Inf),fill="darkred", alpha=0.2) +  geom_rect(data=NULL,aes(xmin=-0.15,xmax=0,ymin=-Inf,ymax=Inf),fill="darkblue",  alpha=0.1)+geom_rect(data=NULL,aes(xmin=-.35,xmax=-.15,ymin=-Inf,ymax=Inf),fill="darkblue",alpha=0.2) + geom_rect(data=NULL,aes(xmin=-0.65,xmax=-.35,ymin=-Inf,ymax=Inf),fill="darkblue",alpha=0.3)  + geom_rect(data=NULL,aes(xmin=-Inf,xmax=-.65,ymin=-Inf,ymax=Inf),fill="darkblue", alpha=0.4) + geom_histogram(aes(x=STIAT.test), data=data, breaks=seq(-1.3, 1.3, by =.05), col="grey30",fill="white", alpha = .75)  + geom_vline(xintercept = -.65, linetype="dashed", size=.3) + geom_vline(xintercept = -.15, linetype="dashed", size=.3) + geom_vline(xintercept = .15, linetype="dashed", size=.3) + geom_vline(xintercept = -.35, linetype="dashed", size=.3) + geom_vline(xintercept = .65, linetype="dashed", size=.3) + geom_vline(xintercept = .35, linetype="dashed", size=.3)+ geom_vline(xintercept = 0) + expand_limits(x=c(-1.15,1.15)) + annotate("text", x = -.95, y = 56, label = "Strong",  size = 2.25) + annotate("text", x = -.95, y = 54, label = "negative",  size = 2.25) + annotate("text", x = .95, y = 56, label = "Strong",  size = 2.25) + annotate("text", x = .95, y = 54, label = "positive",  size = 2.25) + annotate("text", x = -.5, y = 56, label = "Moderate",  size = 2.25) + annotate("text", x = -.5, y =54, label = "negative",  size = 2.25)   + annotate("text", x = .5, y = 56, label = "Moderate",  size = 2.25) + annotate("text", x = .5, y = 54, label = "positive",  size = 2.25)   + annotate("text", x = .25, y = 56, label = "Weak",  size = 2.25) + annotate("text", x = .25, y = 54, label = "positive",  size = 2.25) + annotate("text", x = -.25, y = 56, label = "Weak",  size = 2.25) + annotate("text", x = -.25, y = 54, label = "negative",  size = 2.25) + annotate("text", x = 0, y = 56, label = "Neutral",  size = 2.25) + guides(color=guide_legend(title=NULL)) + theme_bw() + theme(axis.title.x = element_text(size = 12, vjust = .25)) + xlab("Implicit Attitudes Towards Sisi (implicit.sisi)")  + ylab("Count") + guides(color=guide_legend(title=NULL)) + theme_bw() + theme(axis.title.x = element_text(size = 12, vjust = .25))
dev.off()

#Figure 4: Point Estimates

direct.support.est<-mean(data$direct.support,na.rm=TRUE)
direct.support.se<-sd(data$direct.support, na.rm=TRUE)/sqrt(length(data$direct.support[is.na(data$direct.support)==FALSE]))

data$direct.z<- as.numeric(scale(data$approve_abdel_rev, center = TRUE, scale = TRUE))

direct.support.w.est<-weighted.mean(data$direct.support, data$w,na.rm=TRUE)
direct.support.w.se<-sqrt(wtd.var(data$direct.support, data$w,na.rm=TRUE))/sqrt(length(data$direct.support[is.na(data$direct.support)==FALSE]))

implicit.support.est<-mean(data$implicit.support,na.rm=TRUE)
implicit.support.se<-sd(data$implicit.support, na.rm=TRUE)/sqrt(length(data$implicit.support[is.na(data$implicit.support)==FALSE]))

data$implicit.z<- as.numeric(scale(data$STIAT.test, center = TRUE, scale = TRUE))

mean(data$STIAT.test,na.rm=TRUE)
sd(data$STIAT.test,na.rm=TRUE)

implicit.support.w.est<-weighted.mean(data$implicit.support,data$w,na.rm=TRUE)
implicit.support.w.se<-sqrt(wtd.var(data$implicit.support, data$w,na.rm=TRUE))/sqrt(length(data$implicit.support[is.na(data$implicit.support)==FALSE]))

data$list2a.clean[data$list2a.clean==98]<-NA
data$list2b.clean[data$list2b.clean==98]<-NA

mean(data$list2a.clean, na.rm=TRUE) - mean(data$list2b.clean, na.rm=TRUE) 

list.exp2<-data.frame(cbind(data$list2a.clean,data$list2b.clean))
list.exp2$list_2<-rowSums(list.exp2, na.rm=TRUE)
list.exp2<-data.frame(cbind(list.exp2,data$w))

list.exp2<-subset(list.exp2, is.na(X1)==FALSE | is.na(X2)==FALSE)
list.exp2$sisi.treat<-1
list.exp2$sisi.treat[is.na(list.exp2$X2)==FALSE]<-0
list.exp2<-list.exp2[,3:5]
ict.test(list.exp2$list_2, list.exp2$sisi.treat, J = 3, alpha = 0.05, n.draws = 250000, gms = TRUE, pi.table = TRUE)

set.seed(1)
diff.in.means.results <- ictreg(list_2 ~ 1, data = list.exp2, treat = "sisi.treat", J=3, method = "lm")
summary(diff.in.means.results)

set.seed(1)
diff.in.means.results <- ictreg(list_2 ~ 1, data = list.exp2, treat = "sisi.treat", J=3, method = "lm")
summary(diff.in.means.results)
list2.support.est<-summary(diff.in.means.results)$par.treat
list2.support.se<-summary(diff.in.means.results)$se.treat

list.exp2.w<-subset(list.exp2,is.na(list.exp2$data.w)==FALSE)

set.seed(1)
diff.in.means.results <- ictreg(list_2 ~ 1, data = list.exp2.w, treat = "sisi.treat", J=3, method = "lm",weight="data.w")
summary(diff.in.means.results)
list2.support.w.est<-summary(diff.in.means.results)$par.treat
list2.support.w.se<-summary(diff.in.means.results)$se.treat

label<-c("Direct","Direct (wt)","List","List (wt)","Implicit","Implicit (wt)")
mean<-c(direct.support.est,direct.support.w.est,list2.support.est,list2.support.w.est,implicit.support.est,implicit.support.w.est)
sem<-c(direct.support.se,direct.support.w.se,list2.support.se,list2.support.w.se,implicit.support.se,implicit.support.w.se)
estimates<-data.frame(cbind(label, mean, sem)) 
estimates$mean<-as.numeric(paste(estimates$mean)) 
estimates$lci95<-mean-1.96*sem 
estimates$uci95<-mean+1.96*sem 
estimates[4,5]=1.00

pdf('fig4-estimates-plot.pdf', width=6, height=4) 
ggplot(estimates, aes(x=label, y=mean)) + geom_point(size=3, color="darkblue") + geom_errorbar(aes(ymin=lci95, ymax=uci95), width=.1, color="darkblue") + xlab("Measurement Strategy") + ylab("Estimated Proportion of Sisi Supporters") + theme(plot.title = element_text(lineheight=.8, face="bold")) + theme_bw() + theme(legend.title=element_blank()) + geom_hline(yintercept = 0, linetype="dashed", size=.3) + geom_hline(yintercept = 1, linetype="dashed", size=.3) + expand_limits(y=c(0,1)) + scale_y_continuous(limits = c(0, 1.00)) 
dev.off()

#Figure 5: Comparison of D-score and Explicit Measure

data$type<-NA
data$type[data$direct.support==1 & data$implicit.support==0]<-"Classic Attitude Dissociation"
data$type[data$direct.support==0 & data$implicit.support==0]<-"Congruent Negative"
data$type[data$direct.support==1 & data$implicit.support==1]<-"Congruent Positive"
data$type[data$direct.support==0 & data$implicit.support==1]<-"Inverse Attitude Dissociation"
prop.table(ftable(data$type)) 

m1<-zelig(scale(STIAT.test)~scale(approve_abdel_rev), model="ls", data=data) 
summary(m1) 

cor(data$STIAT.test, data$approve_abdel_rev, use="complete.obs")

pdf('fig5-dscorevexplicit-scatter.pdf', width=6, height=6)                                                                                                                                                                                                                                                                                                                                                                                                                                              
ggplot(data, aes(x=STIAT.test, y=approve_abdel_rev)) + xlab("Implicit Attitudes Towards Sisi (implicit.sisi)") + ylab("Explicit Attitudes Towards Sisi (explicit.sisi)") + geom_smooth(method=lm, se=FALSE, colour="black") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE)+ theme_bw() + theme(legend.title=element_blank()) + geom_jitter(color="darkblue", alpha=.25,size=1.25, width=.5, height=1.25) + scale_x_continuous(limits = c(-1, 1))+ scale_y_continuous(limits = c(-.5, 5.5)) + geom_vline(xintercept = 0, linetype="dashed", size=.3) + geom_hline(yintercept = 2.5, linetype="dashed", size=.3) + annotate("text", x = -.5, y = 4, label = "Classic Dissociation",  size = 4.25) + geom_hline(yintercept = 2.5, linetype="dashed", size=.3) + annotate("text", x = .5, y = 1, label = "Inverse Dissociation",  size = 4.25) + annotate("text", x = -.5, y = 3.75, label = "19.1%",  size = 4.25) + annotate("text", x = .5, y = .75, label = "25.3%",  size = 4.25) + annotate("text", x = -.9, y =5, label = "r = 0.17",  size = 4.25)
dev.off()

#Figure 6: Determinants of Expressed Explicit and Implicit Support (standardized approach)
dvlist<-c("rel_christian", "female", "work_govemp", "highinc", "islamist", "lowed", "work_student", "liberal")

estimates <- matrix(data = NA, nrow=96, ncol=8)
estimates<-data.frame(estimates)
colnames(estimates)<-c("iv","dv", "weights", "cov","est","se","t","pvalue")

m1<-lm(direct.z ~ rel_christian, data = data, weights=data$w) 
m1<-lm(direct.z ~ rel_christian,   data = data, weights=data$w)      
m2<-lm(direct.z ~ female,   data = data, weights=data$w)      
m3<-lm(direct.z ~ work_govemp,   data = data, weights=data$w)      
m4<-lm(direct.z ~ highinc,   data = data, weights=data$w)      
m5<-lm(direct.z ~ islamist,   data = data, weights=data$w)      
m6<-lm(direct.z ~ lowed,   data = data, weights=data$w)      
m7<-lm(direct.z ~ work_student,   data = data, weights=data$w)      
m8<-lm(direct.z ~ liberal,   data = data, weights=data$w)      
m9<-lm(implicit.z ~ rel_christian,   data = data, weights=data$w)      
m10<-lm(implicit.z ~ female,   data = data, weights=data$w)      
m11<-lm(implicit.z ~ work_govemp,   data = data, weights=data$w)      
m12<-lm(implicit.z ~ highinc,   data = data, weights=data$w)      
m13<-lm(implicit.z ~ islamist,   data = data, weights=data$w)      
m14<-lm(implicit.z ~ lowed,   data = data, weights=data$w)      
m15<-lm(implicit.z ~ work_student,   data = data, weights=data$w)      
m16<-lm(implicit.z ~ liberal,   data = data, weights=data$w)      
m17<-lm(direct.z ~ rel_christian,   data = data, weights=data$w.tr)      
m18<-lm(direct.z ~ female,   data = data, weights=data$w.tr)      
m19<-lm(direct.z ~ work_govemp,   data = data, weights=data$w.tr)      
m20<-lm(direct.z ~ highinc,   data = data, weights=data$w.tr)      
m21<-lm(direct.z ~ islamist,   data = data, weights=data$w.tr)      
m22<-lm(direct.z ~ lowed,   data = data, weights=data$w.tr)      
m23<-lm(direct.z ~ work_student,   data = data, weights=data$w.tr)      
m24<-lm(direct.z ~ liberal,   data = data, weights=data$w.tr)      
m25<-lm(implicit.z ~ rel_christian,   data = data, weights=data$w.tr)      
m26<-lm(implicit.z ~ female,   data = data, weights=data$w.tr)      
m27<-lm(implicit.z ~ work_govemp,   data = data, weights=data$w.tr)      
m28<-lm(implicit.z ~ highinc,   data = data, weights=data$w.tr)      
m29<-lm(implicit.z ~ islamist,   data = data, weights=data$w.tr)      
m30<-lm(implicit.z ~ lowed,   data = data, weights=data$w.tr)      
m31<-lm(implicit.z ~ work_student,   data = data, weights=data$w.tr)      
m32<-lm(implicit.z ~ liberal,   data = data, weights=data$w.tr)      

for (i in 1:32) {
  eval(parse(text=paste("estimates[",i,",5:8]<- summary(m",i,")$coefficients[2,c(1,2,3,4)]",sep="")))
}

estimates$cov[1:32]<-"bivariate"
estimates$cov[33:64]<-"demographic"
estimates$cov[65:96]<-"region"

estimates$weights[1:16]<-"w"
estimates$weights[17:32]<-"w.tr"
estimates$weights[33:48]<-"w"
estimates$weights[49:64]<-"w.tr"
estimates$weights[65:80]<-"w"
estimates$weights[81:96]<-"w.tr"

estimates$dv[1:8]<-"direct.z"
estimates$dv[9:16]<-"implicit.z"
estimates$dv[17:24]<-"direct.z"
estimates$dv[25:32]<-"implicit.z"
estimates$dv[33:40]<-"direct.z"
estimates$dv[41:48]<-"implicit.z"
estimates$dv[49:56]<-"direct.z"
estimates$dv[57:64]<-"implicit.z"
estimates$dv[65:72]<-"direct.z"
estimates$dv[73:80]<-"implicit.z"
estimates$dv[81:88]<-"direct.z"
estimates$dv[89:96]<-"implicit.z"

estimates$iv<-rep(dvlist, times=12)

m33<-lm(direct.z ~ rel_christian + female + work_govemp + highinc + islamist + lowed + work_student + liberal + age,   data = data, weights=data$w)      
estimates[33,5:8]<- summary(m33)$coefficients[2,c(1,2,3,4)]
estimates[34,5:8]<- summary(m33)$coefficients[3,c(1,2,3,4)]
estimates[35,5:8]<- summary(m33)$coefficients[4,c(1,2,3,4)]
estimates[36,5:8]<- summary(m33)$coefficients[5,c(1,2,3,4)]
estimates[37,5:8]<- summary(m33)$coefficients[6,c(1,2,3,4)]
estimates[38,5:8]<- summary(m33)$coefficients[7,c(1,2,3,4)]
estimates[39,5:8]<- summary(m33)$coefficients[8,c(1,2,3,4)]
estimates[40,5:8]<- summary(m33)$coefficients[9,c(1,2,3,4)]

m34<-lm(implicit.z ~ rel_christian + female + work_govemp + highinc + islamist + lowed + work_student + liberal + age,   data = data, weights=data$w)      
estimates[41,5:8]<- summary(m34)$coefficients[2,c(1,2,3,4)]
estimates[42,5:8]<- summary(m34)$coefficients[3,c(1,2,3,4)]
estimates[43,5:8]<- summary(m34)$coefficients[4,c(1,2,3,4)]
estimates[44,5:8]<- summary(m34)$coefficients[5,c(1,2,3,4)]
estimates[45,5:8]<- summary(m34)$coefficients[6,c(1,2,3,4)]
estimates[46,5:8]<- summary(m34)$coefficients[7,c(1,2,3,4)]
estimates[47,5:8]<- summary(m34)$coefficients[8,c(1,2,3,4)]
estimates[48,5:8]<- summary(m34)$coefficients[9,c(1,2,3,4)]

m35<-lm(direct.z ~ rel_christian + female + work_govemp + highinc + islamist + lowed + work_student + liberal + age,   data = data, weights=data$w.tr)      
estimates[49,5:8]<- summary(m35)$coefficients[2,c(1,2,3,4)]
estimates[50,5:8]<- summary(m35)$coefficients[3,c(1,2,3,4)]
estimates[51,5:8]<- summary(m35)$coefficients[4,c(1,2,3,4)]
estimates[52,5:8]<- summary(m35)$coefficients[5,c(1,2,3,4)]
estimates[53,5:8]<- summary(m35)$coefficients[6,c(1,2,3,4)]
estimates[54,5:8]<- summary(m35)$coefficients[7,c(1,2,3,4)]
estimates[55,5:8]<- summary(m35)$coefficients[8,c(1,2,3,4)]
estimates[56,5:8]<- summary(m35)$coefficients[9,c(1,2,3,4)]

m36<-lm(implicit.z ~ rel_christian + female + work_govemp + highinc + islamist + lowed + work_student + liberal + age,   data = data, weights=data$w.tr)      
estimates[57,5:8]<- summary(m36)$coefficients[2,c(1,2,3,4)]
estimates[58,5:8]<- summary(m36)$coefficients[3,c(1,2,3,4)]
estimates[59,5:8]<- summary(m36)$coefficients[4,c(1,2,3,4)]
estimates[60,5:8]<- summary(m36)$coefficients[5,c(1,2,3,4)]
estimates[61,5:8]<- summary(m36)$coefficients[6,c(1,2,3,4)]
estimates[62,5:8]<- summary(m36)$coefficients[7,c(1,2,3,4)]
estimates[63,5:8]<- summary(m36)$coefficients[8,c(1,2,3,4)]
estimates[64,5:8]<- summary(m36)$coefficients[9,c(1,2,3,4)]

m37<-lm(direct.z ~ rel_christian + female + work_govemp + highinc + islamist + lowed + work_student + liberal + age + region.indupper + region.inddelta + region.indalex + region.inddesert,   data = data, weights=data$w)      
estimates[65,5:8]<- summary(m37)$coefficients[2,c(1,2,3,4)]
estimates[66,5:8]<- summary(m37)$coefficients[3,c(1,2,3,4)]
estimates[67,5:8]<- summary(m37)$coefficients[4,c(1,2,3,4)]
estimates[68,5:8]<- summary(m37)$coefficients[5,c(1,2,3,4)]
estimates[69,5:8]<- summary(m37)$coefficients[6,c(1,2,3,4)]
estimates[70,5:8]<- summary(m37)$coefficients[7,c(1,2,3,4)]
estimates[71,5:8]<- summary(m37)$coefficients[8,c(1,2,3,4)]
estimates[72,5:8]<- summary(m37)$coefficients[9,c(1,2,3,4)]

m38<-lm(implicit.z ~ rel_christian + female + work_govemp + highinc + islamist + lowed + work_student + liberal + age + region.indupper + region.inddelta + region.indalex + region.inddesert,   data = data, weights=data$w)      
estimates[73,5:8]<- summary(m38)$coefficients[2,c(1,2,3,4)]
estimates[74,5:8]<- summary(m38)$coefficients[3,c(1,2,3,4)]
estimates[75,5:8]<- summary(m38)$coefficients[4,c(1,2,3,4)]
estimates[76,5:8]<- summary(m38)$coefficients[5,c(1,2,3,4)]
estimates[77,5:8]<- summary(m38)$coefficients[6,c(1,2,3,4)]
estimates[78,5:8]<- summary(m38)$coefficients[7,c(1,2,3,4)]
estimates[79,5:8]<- summary(m38)$coefficients[8,c(1,2,3,4)]
estimates[80,5:8]<- summary(m38)$coefficients[9,c(1,2,3,4)]

m39<-lm(direct.z ~ rel_christian + female + work_govemp + highinc + islamist + lowed + work_student + liberal + age + region.indupper + region.inddelta + region.indalex + region.inddesert,   data = data, weights=data$w.tr)      
estimates[81,5:8]<- summary(m39)$coefficients[2,c(1,2,3,4)]
estimates[82,5:8]<- summary(m39)$coefficients[3,c(1,2,3,4)]
estimates[83,5:8]<- summary(m39)$coefficients[4,c(1,2,3,4)]
estimates[84,5:8]<- summary(m39)$coefficients[5,c(1,2,3,4)]
estimates[85,5:8]<- summary(m39)$coefficients[6,c(1,2,3,4)]
estimates[86,5:8]<- summary(m39)$coefficients[7,c(1,2,3,4)]
estimates[87,5:8]<- summary(m39)$coefficients[8,c(1,2,3,4)]
estimates[88,5:8]<- summary(m39)$coefficients[9,c(1,2,3,4)]

m40<-lm(implicit.z ~ rel_christian + female + work_govemp + highinc + islamist + lowed + work_student + liberal + age + region.indupper + region.inddelta + region.indalex + region.inddesert,  data = data, weights=data$w.tr)      
estimates[89,5:8]<- summary(m40)$coefficients[2,c(1,2,3,4)]
estimates[90,5:8]<- summary(m40)$coefficients[3,c(1,2,3,4)]
estimates[91,5:8]<- summary(m40)$coefficients[4,c(1,2,3,4)]
estimates[92,5:8]<- summary(m40)$coefficients[5,c(1,2,3,4)]
estimates[93,5:8]<- summary(m40)$coefficients[6,c(1,2,3,4)]
estimates[94,5:8]<- summary(m40)$coefficients[7,c(1,2,3,4)]
estimates[95,5:8]<- summary(m40)$coefficients[8,c(1,2,3,4)]
estimates[96,5:8]<- summary(m40)$coefficients[9,c(1,2,3,4)]

estimates$l95ci<-estimates$est-1.96*estimates$se
estimates$u95ci<-estimates$est+1.96*estimates$se

estimates$group<-NA
estimates$group[estimates$dv=="direct.z" & estimates$cov=="bivariate"]<-1
estimates$group[estimates$dv=="direct.z" & estimates$cov=="demographic"]<-2
estimates$group[estimates$dv=="direct.z" & estimates$cov=="region"]<-3
estimates$group[estimates$dv=="implicit.z" & estimates$cov=="bivariate"]<-4
estimates$group[estimates$dv=="implicit.z" & estimates$cov=="demographic"]<-5
estimates$group[estimates$dv=="implicit.z" & estimates$cov=="region"]<-6

estimates$varname[estimates$iv=="female"]<-"Female"
estimates$varname[estimates$iv=="work_student"]<-"Student"
estimates$varname[estimates$iv=="lowed"]<-"Low Education"
estimates$varname[estimates$iv=="highinc"]<-"High Income"
estimates$varname[estimates$iv=="work_govemp"]<-"Govt Emp"
estimates$varname[estimates$iv=="rel_christian"]<-"Christian"
estimates$varname[estimates$iv=="islamist"]<-"Islamist"
estimates$varname[estimates$iv=="liberal"]<-"Liberal"
write.csv(estimates,"estimates.csv")

estimates.w<-subset(estimates, estimates$weights=="w")
estimates.w.tr<-subset(estimates, estimates$weights=="w.tr")

pdf('fig6-effects-w.pdf', width=9, height=4.5) 
ggplot(estimates.w, aes(x=varname, y=est, group=group,colour=dv, shape=cov)) + geom_point(size=3,position=position_dodge(width=.5)) + geom_errorbar(aes(ymin=l95ci, ymax=u95ci), position=position_dodge(width=.5),width=0.2) + xlab("Independent Variable") + ylab("Coefficient Estimate") + theme(plot.title = element_text(lineheight=.8, face="bold")) + theme_bw()  + geom_hline(yintercept = 0, linetype="dashed", size=.3) + scale_y_continuous(limits = c(-1.2, 1.2)) + scale_colour_manual(name="Dependent Variable", values=c("darkred", "darkblue"),labels=c("direct.sisi.z", "implicit.sisi.z")) + theme(legend.position = "right") + scale_shape_discrete(name="Covariates",labels=c("Bivariate", "Demographic","Demographic + Region"))
dev.off() 

pdf('fig6-effects-wtr.pdf', width=9, height=4.5)  
ggplot(estimates.w.tr, aes(x=varname, y=est, group=group,colour=dv, shape=cov)) + geom_point(size=3,position=position_dodge(width=.5)) + geom_errorbar(aes(ymin=l95ci, ymax=u95ci), position=position_dodge(width=.5),width=0.2) + xlab("Independent Variable") + ylab("Coefficient Estimate") + theme(plot.title = element_text(lineheight=.8, face="bold")) + theme_bw()  + geom_hline(yintercept = 0, linetype="dashed", size=.3) + scale_y_continuous(limits = c(-1.5, 1.5)) + scale_colour_manual(name="Dependent Variable", values=c("darkred", "darkblue"),labels=c("direct.sisi.z", "implicit.sisi.z")) + theme(legend.position = "right") + scale_shape_discrete(name="Covariates",labels=c("Bivariate", "Demographic","Demographic + Region"))
dev.off() 

##SUMMARY STATISTICS FOR IAT##
#Error rates for Implicit Association Test

included<-data.frame(c(data$session_id))
included$included<-1
colnames(included)<-c("session_id","included")

iat<-merge(iat,included,by="session_id",all=TRUE)
iat.included<-subset(iat, iat$included==1)
iat.included<-subset(iat.included, iat.included$block_number!=0 & is.na(iat.included$block_number)==FALSE)

table(iat$session_id,iat$block_number)

iat.included$trial_error<-as.numeric(paste(iat.included$trial_error))
iat.included$trial_latency<-as.numeric(paste(iat.included$trial_latency))

iat.included$latency10000<-0
iat.included$latency10000[iat.included$trial_latency>10000]<-1

iat.included$latency400<-0
iat.included$latency400[iat.included$trial_latency<400]<-1

mean(iat.included$latency10000, na.rm=TRUE)
mean(iat.included$latency400, na.rm=TRUE)


#Figure A4: Distribution of Response Times#
pdf('figa4-responsetimes-density.pdf', width=6, height=4)                                                                                                                                                                                                                                                                                                                                                                                                                                              
ggplot() + geom_line(aes(x=trial_latency), trim=TRUE, adjust=2.5, size=1, stat="density", data=iat.included)  + guides(color=guide_legend(title=NULL)) + theme_bw() + theme(axis.title.x = element_text(size = 12, vjust = .25)) + xlab("Reaction Time (milliseconds)") + ylab("Density") + geom_vline(xintercept = 10000, linetype="dashed", size=.3) + geom_vline(xintercept = 400, linetype="dashed", size=.3) + scale_x_continuous(limits = c(0, 12000))
dev.off()

#Figure A5: Distribution of IAT Completion Times#
iat.stats<-ddply(iat, .(session_id), summarize,  Rate1=mean(trial_error), Rate2=mean(trial_latency), Rate3=sum(trial_latency))
colnames(iat.stats)<-c("session_id","mean.trial_error","mean.trial_latency","total.trial_latency")
iat.stats$total.iat.secs<-iat.stats$total.trial_latency/1000 
iat.stats$total.iat.mins<-iat.stats$total.iat.secs/60
mean(iat.stats$total.iat.mins)

iat.stats<-merge(iat.stats,included,by="session_id",all=TRUE)
iat.stats.included<-subset(iat.stats, iat.stats$included==1)

pdf('figa5-iatcompletiontimes-density.pdf', width=6, height=4)                                                                                                                                                                                                                                                                                                                                                                                                                                              
ggplot() + geom_line(aes(x=total.iat.mins), trim=TRUE, adjust=2.5, size=1, stat="density", data=iat.stats.included)  + guides(color=guide_legend(title=NULL)) + theme_bw() + theme(axis.title.x = element_text(size = 12, vjust = .25)) + xlab("Total Time to Complete IAT (minutes)") + ylab("Density") + scale_x_continuous(limits = c(0, 15))
dev.off()

#Figure A6: Error Rates and Response Times Across Blocks#

mean(iat.included$trial_error, na.rm=TRUE) 
mean(iat.included$trial_error[iat.included$block_number==1], na.rm=TRUE)
mean(iat.included$trial_error[iat.included$block_number==2], na.rm=TRUE)
mean(iat.included$trial_error[iat.included$block_number==3], na.rm=TRUE)
mean(iat.included$trial_error[iat.included$block_number==4], na.rm=TRUE)
mean(iat.included$trial_error[iat.included$block_number==5], na.rm=TRUE)

mean(iat.included$trial_latency, na.rm=TRUE) 
mean(iat.included$trial_latency[iat.included$block_number==1], na.rm=TRUE)
mean(iat.included$trial_latency[iat.included$block_number==2], na.rm=TRUE)
mean(iat.included$trial_latency[iat.included$block_number==3], na.rm=TRUE)
mean(iat.included$trial_latency[iat.included$block_number==4], na.rm=TRUE)
mean(iat.included$trial_latency[iat.included$block_number==5], na.rm=TRUE)

blocks <- ddply(iat.included, "block_number", summarise, trial_latency = mean(trial_latency), error=mean(trial_error))
blocks$type<-c("practice","practice","test","practice","test")

pdf('figa6-errorrates-barplot.pdf', width=6, height=4)  
plot1<-ggplot(blocks, aes(x = factor(block_number), y = trial_latency)) + geom_bar(stat = "identity", colour="darkblue",fill="darkblue", alpha=.5) + theme_bw() + theme(axis.title.x = element_text(size = 12, vjust = .25)) + xlab("Block Number") + ylab("Mean Reaction Time (milliseconds)")
plot2<-ggplot(blocks, aes(x = factor(block_number), y = error)) + geom_bar(stat = "identity", colour="darkblue",fill="darkblue", alpha=.5) + theme_bw() + theme(axis.title.x = element_text(size = 12, vjust = .25)) + xlab("Block Number") + ylab("Mean Error Rate")
arrange_ggplot2(plot1,plot2,nrow=1,ncol=2)
dev.off()







  


